/*

  Simple cross-product, outputs in correct triangular form.

  - Coalescing memory access in the Column reads
  - No memory coalescing in the Row reads
  - 50% coalescing in the write (float8s)
  - All threads read in all data

 */

CUBE_KERNEL(simple1x1, cuFloatComplex *array_d, cuFloatComplex *product_d, 
	    const int write)
{
  CUBE_START;

  //set frequency number from blockIdx.y
  int f = blockIdx.y;

  int Row, Col, blockX, blockY;
  CUBE_DEVICE_CALL(findPosition, Col, Row, blockX, blockY);

  //instantiate sum variables
  float sumXXreal = 0.0, sumXXimag = 0.0;
  float sumXYreal = 0.0, sumXYimag = 0.0;
  float sumYXreal = 0.0, sumYXimag = 0.0;
  float sumYYreal = 0.0, sumYYimag = 0.0;

  for(int t=0; t<NTIME; t++){
    float4 col = ((float4*)array_d)[t*NFREQUENCY*NSTATION + f*NSTATION + Col];
    float4 row = ((float4*)array_d)[t*NFREQUENCY*NSTATION + f*NSTATION + Row];
    CUBE_ADD_BYTES(32);

    sumXXreal += row.x * col.x;
    sumXXreal += row.y * col.y;
    sumXXimag += row.y * col.x;
    sumXXimag -= row.x * col.y;
    sumXYreal += row.x * col.z;
    sumXYreal += row.y * col.w;
    sumXYimag += row.y * col.z;
    sumXYimag -= row.x * col.w;
    sumYXreal += row.z * col.x;
    sumYXreal += row.w * col.y;
    sumYXimag += row.w * col.x;
    sumYXimag -= row.z * col.y;
    sumYYreal += row.z * col.z;
    sumYYreal += row.w * col.w;
    sumYYimag += row.w * col.z;
    sumYYimag -= row.z * col.w;
  } 

  if (write)
    CUBE_DEVICE_CALL(write1x1, Col, Row, product_d,  
		     sumXXreal,  sumXXimag, sumXYreal,  sumXYimag,  
		     sumYXreal, sumYXimag,  sumYYreal,  sumYYimag);

  if (Col <= Row) {
    CUBE_ADD_BYTES(32);         //output bytes
    CUBE_ADD_FLOPS(NTIME*32);
  }

  CUBE_END;
}
